
************************************************
**** Do file to create figures in paper
**** The life, death, and diversity of pro-government militias:
**** The fully revised pro-government militias database version 2.0
**** Sabine C. Carey, Neil J. Mitchell and Katrin Paula
**** Research & Politics
**** October 2021
**** Stata version: 17.0
************************************************


********************************************************************************
*
* Packages: Before running this code, please make sure you have the following
*			.ado-files installed
*
*			- grc1leg Package (user written)
*					net describe grc1leg, from(http://www.stata.com/users/vwiggins)
*					net install grc1leg.pkg
*
*			- spineplot Package
*					ssc install spineplot
*
*
********************************************************************************
********************************************************************************
* Content:	Replication of results in manuscript and supporting information	
*
*			Steps:
*			1) Set working directory
*			2) Run Do-File, it will
*				a) Automatically load the necessary dataset
*				b) Replicate results and figures of manuscript
*			
*
********************************************************************************


	set more off
	clear all
  
********************************************************************************
* Setting working directory 
********************************************************************************
	
	cd "/replication_material/" //Set path for your working directory here
		* Example for Mac: cd "/Users/name/Dektop/Replication" 
		* Example for Windows "C:\Users\name\Desktop\Replication" 
	

********************************************************************************
* Producing log-file of results
********************************************************************************

	* starting log file
	log using "CareyMitchellPaula_replication_RP", replace 
	

***********************************************
* Figures using country-year format 
***********************************************

	* load data in country-year format
	use "pgmdv2_countryyear.dta"
	
	* aggregate data to year-format
	
		* all pgms per year
		sort year
		bysort year: egen presence_sum=total(presence_count)
	
		* informal pgms per year
		bysort  year: egen presence_informal_sum=total(presence_informal_count) 
		
		* semi-officials pgms per year
		bysort  year: egen presence_semi_sum=total(presence_semiofficial_count) 
		
		* countries with PGM presence per year
		bysort year: egen presence_country = total(presence)
		
		* countries with PGMs fromed by government per year
		bysort year: egen presence_gov_formed_country = total(presence_govformed)
		
		* countries with newly government-formed PGMs per year
		bysort year: egen presence_gov_formed_on_country = total(presence_gov_formed_on)
		
	* collapse data for plotting
	collapse (first) presence_sum presence_informal_sum presence_semi_sum presence_country presence_gov_formed_country presence_gov_formed_on_country, by(year)
	
	

* Figure 1. Number of PGMs over time
*--------------------------------------------------

	#delimit ;
		twoway line presence_sum presence_informal_sum presence_semi_sum year,
		sort 
		lcolor(black ebblue emidblue)
		xline(0, lstyle(grid) lcolor(white) lwidth(*1.5))
		xlabel(1981 1985(5)2010 2014, grid glcolor(white)) 
		ylabel(,nogrid)
		xtitle("")
		ytitle("Number of PGMs")
		xline(2007, lpattern(shortdash) lcolor(black)) 
		legend(order (1 "Any PGM" 2 "Informal PGMs" 3 "Semi-official PGMs") 
		size(small) keygap(*0.5) symysize(*0.5) region(lcolor(white)) rows(1) cols(3))
		note("Excluding Somalia 1981-2011 & Lebanon 1981-2007")
		graphregion(color(white))
	;
	#delimit cr  	
	graph export "figures/Figure1.pdf", replace
	

* Figure 2. Distribution of PGMs across countries and time
*--------------------------------------------------

	#delimit ;
		twoway line presence_country presence_gov_formed_country presence_gov_formed_on_country year, 
		sort 
		xlabel(1981 1985(5)2010 2014) 	
		title("Countries with PGMs")
		ytitle("Number of countries with PGMs")
		title("")
		ylabel(0(10)60,nogrid)
		xline(2007, lpattern(shortdash) lcolor(black)) 
		lcolor(black ebblue eltgreen)
		legend(order (1 "Any PGM" 2 "PGM formed by government" 3 "New PGM formed by government") 
		size(small) keygap(*0.5) symysize(*0.5) region(lcolor(white)) rows(3) cols(1)) 	
		note("Excluding Somalia 1981-2011 & Lebanon 1981-2007")
		graphregion(color(white))
	;
	#delimit cr  
	graph export "figures/Figure2.pdf", replace 
	


***********************************************
* Figures using group format 
***********************************************
	
	* load data in group format
	clear all
	use "pgmdv2_group.dta"

	
	
* Figure 3. Primary characteristic of PGM membership, by government formation
*--------------------------------------------------
			
	* create dummies for primary membership plotting
	
		* for ethnicity
		gen prim_ethnicity = primary_mem
		recode prim_ethnicity (0=.)
		recode prim_ethnicity (9=.)
		recode prim_ethnicity (2/8=0) 
		recode prim_ethnicity (0=1) if alt_primary_mem==1 
	
		* for local
		gen prim_local = primary_mem
		recode prim_local (0=.)
		recode prim_local (9=.)
		recode prim_local (1/2=0) (3=1) (4/8=0)
		recode prim_local (0=1) if alt_primary_mem==3
	
		* for non-civilian
		gen prim_non_civilian = primary_mem
		recode prim_non_civilian  (0=.)
		recode prim_non_civilian  (9=.)
		recode prim_non_civilian (1/4=0) (5=1) (6/8=0)
		recode prim_non_civilian (0=1) if alt_primary_mem==5
	
		* for political
		gen prim_political = primary_mem
		recode prim_political (0=.)
		recode prim_political (9=.)
		recode prim_political (1/5=0) (6=1) (7/8=0)
		recode prim_political (0=1) if alt_primary_mem==6
	
		* for religious
		gen prim_reli = primary_mem
		recode prim_reli (0=.)
		recode prim_reli (9=.)
		recode prim_reli (1/6=0) (7=1) (8=0)
		recode prim_reli (0=1) if alt_primary_mem==7
	
		* for other (incl. other/nationalist/ideological)
		gen prim_other = primary_mem
		recode prim_other (0=.)
		recode prim_other (9=.)
		recode prim_other (1=0) (2=1) (3=0) (4=1) (5/7=0) (8=1)
		recode prim_other (0=1) if alt_primary_mem==2 | alt_primary_mem==4
	
	
	* distribution pgms used in graph
	gen diff= prim_ethnicity + prim_local + prim_non_civilian  + prim_political + prim_reli + prim_other
	tab diff gov_formed if gov_formed <3 // tab cases with clear government formation (gov_formed < 3)
		* 195 not formed by government
		* 154 formed by government
		* 38 with 2 primary memberships
		
	* get percentages 
	tab prim_ethnicity gov_formed if gov_formed <3, col
	tab prim_local gov_formed if gov_formed <3, col
	tab prim_non_civilian gov_formed if gov_formed <3, col
	tab prim_reli gov_formed if gov_formed <3, col
	tab prim_political gov_formed if gov_formed <3, col
	tab prim_other gov_formed if gov_formed <3, col
	
		* not formed by government: 
			* 30.26% ethnic, 15.38% local, 14.87% non-civilian, 20.00% religious, 16.41% political, 16.92% other
		* formed by government: 
			* 12.99% ethnic; 37.01% local, 16.23% non-civilian, 7.14% religious, 13.64% political, 20.13% other

	* Figure 3
	#delimit ;
		graph hbar prim_ethnicity prim_local prim_non_civilian  prim_reli prim_political prim_other if gov_formed <3 , 
		over(gov_formed, relabel(1 "No" 2 "Yes")) stack percent
		graphregion(color(white))
		ylabel( none, nogrid)
		yscale(lstyle(none))
		ytitle("")
		bar(1, color(black*0.95)) bar(2, color(black*.73)) bar(3, color(black*0.55)) bar(4, color(black*0.40)) bar(5, color(black*0.3))
		bar(6, color(black*0.25))
		legend(row(2) col(3))
		legend(lab(1 "Ethnic") lab(2 "Local") lab(3 "Non-civilian") lab(4 "Religious") lab(5 "Political") lab(6 "Other")) 
		legend(size(small) position(7) bmargin(small)  region(lcolor(white)))
		ttext(3 23   "12.99%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(27 23  "37.01%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(51 23  "16.23%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(63 23   "7.14%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(72 23  "13.64%" "", place(e) size(vsmall) orientation(horizontal) col(white)) 
		ttext(88 23  "20.13%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(11 78  "30.26%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(30 78  "15.38%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(44 78  "14.87%" "", place(e) size(vsmall) orientation(horizontal) col(white))  
		ttext(59 78  "20.00%" "", place(e) size(vsmall) orientation(horizontal) col(white))
		ttext(75 78  "16.41%" "", place(e) size(vsmall) orientation(horizontal) col(white)) 
		ttext(89 78  "16.92%" "", place(e) size(vsmall) orientation(horizontal) col(white)) 
	;
	# delimit cr
	graph export "figures/Figure3.pdf", replace 

	
	
* Figure 4. Types of violence committed, three most common combinations
*--------------------------------------------------
	
	* get most combinations 
		
		* informal PGMs
		egen group1 = group(violencetype_kill violencetype_kidnap violencetype_torture violencetype_sexual violencetype_beat) if government_relation==1, label
		tab group1 if group1 !=1 //only pgms with any reported violence type at all
		
		* most common three
			* 1 killing: 33.82%  						(10000)
			* 2 killing & kidnapping: 8.21% 			(11000)
			* 3 killing, kidnapping, torture: 8.21%		(11100)
		
		
		* semi-official PGMs
		egen group2 = group(violencetype_kill violencetype_kidnap violencetype_torture violencetype_sexual violencetype_beat) if government_relation==2, label
		tab group2 if group2!=1 //only pgms with any reported violence type at all
		
		* most common three
			* 1 killing: 32.56%  						(10000)
			* 2 killing & beating: 11.36% 				(10001)
			* 3 beating: 8.14%							(00001)
	
	* pre-step 1: create indicator variable for most common combinations
	gen cat1 = 0
	replace cat1 = 1 if group2== 7  // most common combi semi-officials: killing
	replace cat1 = 2 if group2== 8  // 2nd most common semi-officials: killing & beating
	replace cat1 = 3 if group2== 2  // 3rd most common semi-officials: beating
	replace cat1 = 1 if group1== 9  // most common combi informals: killing
	replace cat1 = 2 if group1== 17 // 2nd most common combi informals: killing & kidnapping
	replace cat1 = 3 if group1== 21 // 3rd most common combi informals: killing, kidnapping & torture
	replace cat1=. if group1==. & group2==.
	replace cat1 =. if group1 ==1 | group2==1 // excluding pgms with no reported form of violence
	
	* pre-step 2: create dummy variables for plotting
	tab cat1, gen(dum)
	
	* plot and edit graph to hide dummy1 bar = others (necessary, otherwise % wrong)
	#delimit ;
		graph hbar dum2 dum3 dum4 dum1  if government_relation !=0 , 
		over(government_relation, relabel(1 "Informal" 2 "Semi-official"))  percent
		graphregion(color(white))
		ylabel( none, nogrid)
		yscale(lstyle(none))
		ytitle("")
		title ("")
		bar(1, color(black*0.9)) bar(2, color(black*.6)) bar(3, color(black*0.4)) bar(4, color(black*0.28))
		legend(off)
		ttext(1 20  "8.14%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 31  "11.63%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 41  "32.56%" "", place(e) size(small) orientation(horizontal) col(white) ) 
		
		ttext(35 20  "Beating", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(35 31  "Killing & Beating", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(35 41  "Killing", place(e) size(small) orientation(horizontal) col(black) ) 
		
		ttext(1 71  "8.21%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 82  "8.21%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 93  "33.82%" "", place(e) size(small) orientation(horizontal) col(white) )  
		
		ttext(35 71  "Killing, Kidnapping & Torture", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(35 82  "Killing & Kidnapping", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(35 93  "Killing", place(e) size(small) orientation(horizontal) col(black) )  
	;
	# delimit cr
	
	gr_edit .plotregion1.bars[5].draw_view.setstyle, style(no)
	gr_edit .plotregion1.bars[1].draw_view.setstyle, style(no)
	graph export  "figures/Figure4.pdf", replace
	
		
		
* Figure 5. Reported violence committed by PGM with or without forced membership
*--------------------------------------------------	 
	 
	 * get percentages
	 tab  violencetype_sexual members_coerced if members_coerced < 3, col
	 tab  violencetype_kill members_coerced if members_coerced < 3, col
	 tab  violencetype_beat members_coerced if members_coerced < 3, col
	
	* label for plotting
	label define coerced 1 "Without coercion" 2 "With coercion"
	label value members_coerced coerced
	
	*  Figure 5a. sexual violence
	#delimit ;
		spineplot  violencetype_sexual members_coerced if members_coerced <3,
		bar1(color(black*0.9)) bar2(color(black*0.6)) percent
		graphregion(color(white))
		title("Sexual Violence" "" "", size(large) color(black))
		xti(% Groups using coercion, axis(1) size(medlarge)) 
		xla(0(25)100, axis(1) labsize(medlarge)) 
		yti(% Violence Type Sexual Violence, axis(2) size(medlarge)) 
		yla(0(25)100, axis(2) labsize(medlarge)) 
		legend(off)
		ttext(35 40  "89.86%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(35 88  "62.71%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(95 40  "10.14%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(95 88  "37.29%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		xti("", axis(2) size(medlarge)) 
		xla(, axis(2) labsize(medlarge)) 
       
	;
	# delimit cr
	graph save Figure5a.gph, replace
	
	*  Figure 5b. killings
	#delimit ;
		spineplot  violencetype_kill members_coerced if members_coerced <3,
		bar1(color(black*0.9)) bar2(color(black*0.6)) percent
		graphregion(color(white))
		title("Killings" "" "", size(large) color(black))
		xti(% Groups using coercion, axis(1) size(medlarge)) 
		xla(0(25)100, axis(1) labsize(medlarge)) 
		yti(% Violence Type Killing, axis(2) size(medlarge)) 
		yla(0(25)100, axis(2) labsize(medlarge)) 
		legend(off)
		ttext(25 40  "50.69%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(25 88  "38.98%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(75 40  "49.31%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(75 88  "61.02%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		xti("", axis(2) size(medlarge)) 
		xla(, axis(2) labsize(medlarge)) 
       
	;
	# delimit cr		
	graph save Figure5b.gph, replace
	
	* Figure 5c. beatings
	
		* label variable newly for legend in combined plot
		label def combi 1 "Violence type used" 0"Violence type not used"
		label val violencetype_beat combi
	
	#delimit ;
		spineplot  violencetype_beat members_coerced if members_coerced <3,
		bar1(color(black*0.9)) bar2(color(black*0.6)) percent
		graphregion(color(white))
		title("Beatings" "" "", size(large) color(black))
		xti(% Groups using coercion, axis(1) size(medlarge)) 
		xla(0(25)100, axis(1) labsize(medlarge)) 
		yti(% Violence Type Beatings, axis(2) size(medlarge)) 
		yla(0(25)100, axis(2) labsize(medlarge)) 
		legend(size(small) bmargin(small)  region(lcolor(white)))
		ttext(40 40  "81.80%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(40 88  "69.49%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(90 40  "18.20%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		ttext(90 88  "30.51%" "", place(e) size(medlarge) orientation(horizontal) col(white) )
		xti("", axis(2) size(medlarge)) 
		xla(, axis(2) labsize(medlarge)) 
       
	;
	# delimit cr
	graph save Graph Figure5c.gph, replace
	

	* combine 3 figures into panel
	#delimit ;
		grc1leg Figure5a.gph Figure5b.gph Figure5c.gph, 
		ycommon xcommon row(1)
		graphregion(color(white)) legendfrom(Figure5c.gph)
	;
	# delimit cr
	
	graph display, ysize(5) xsize(17)
	graph export  "figures/Figure5.pdf", replace
	
	erase "Figure5a.gph" // clean-up: erase temporary files from working directory
	erase "Figure5b.gph" 
	erase "Figure5c.gph"
			
			
* Figure 6. Pedigree of PGMs
*--------------------------------------------------	 
	
	* get percentages 
	tab prior_group, m // 61.7; 2.4; 35.9
	tab prior_armed if prior_group==2, m // 10.5; 8.8; 80.7
	tab prior_rebel if prior_armed==2, m // 54.8; 4.1; 41.1
	
	* Figure 6 done by hand

	
* Figure 7. Three most common causes of PGM termination, by government link 
*--------------------------------------------------	

	* generate indicator for PGMs not terminated yet
	gen termtype_not_terminated = 0
	replace termtype_not_terminated = 1 if year_terminated ==.  
		
	* get combinations
	
		* informal PGMs
		egen term1 = group(termtype_disarmed termtype_gov_change termtype_gov_defect termtype_border_change termtype_pgm_defect termtype_integrate termtype_unclear) if termtype_not_terminated ==0 & government_relation ==1, label
		tab term1
	
		* most common three
			* 1 government changed: 45.13%  (0100000)
			* 2 government defected: 19.49% (0010000)
			* 3 border change: 8.21%		(0001000)
		
		* semi-official PGMs
		egen term2 = group(termtype_disarmed termtype_gov_change termtype_gov_defect termtype_border_change termtype_pgm_defect termtype_integrate termtype_unclear) if termtype_not_terminated ==0 & government_relation ==2, label
		tab term2
	
		* most common three
			* 1 government changed: 31.40% 	(0100000)
			* 2 government defected: 26.74% (0010000)
			* 3 PGM integrated: 10.47 %		(0000010)

		
	* pre-step 1: create indicator variable 
	gen kat1 = 0
	replace kat1 =1 if term2== 8  // most common semi-official: government changed
	replace kat1 =2 if term2== 5  // 2nd most common semi-official: government defected
	replace kat1 =3 if term2== 2  // 3rd most common semi-official: PGM integrated
	replace kat1 =1 if term1== 10 // most common informal: government changed
	replace kat1 =2 if term1== 7  // 2nd most common informal: government defected
	replace kat1 =3 if term1== 6  // 3rd most common informal: Border changed
	replace kat1=. if term2==. & term1==.
	
	* pre-Step 2: create dummy variables for plotting
	tab kat1, gen(tum)
	

	* plot and edit graph to hide dummy1 bar = others (necessary, otherwise % wrong)
	#delimit ;
		graph hbar tum2 tum3 tum4 tum1  if government_relation !=0 , 
		over(government_relation, relabel(1 "Informal" 2 "Semi-official"))  percent
		graphregion(color(white))
		ylabel( none, nogrid)
		yscale(lstyle(none))
		ytitle("")
		bar(1, color(black*0.9)) bar(2, color(black*.6)) bar(3, color(blackue*0.4)) bar(4, color(blackue*0.28))
		legend(off)
		ttext(1 20  "10.47%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 31  "26.74%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 41  "31.40%" "", place(e) size(small) orientation(horizontal) col(white) ) 
		
		ttext(48 20  "PGM integreated", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(48 31  "Government defected", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(48 41  "Government changed", place(e) size(small) orientation(horizontal) col(black) ) 
		
		ttext(1 71  "8.21%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 82  "19.49%" "", place(e) size(small) orientation(horizontal) col(white) )  
		ttext(1 93  "45.13%" "", place(e) size(small) orientation(horizontal) col(white) )  
		
		ttext(48 71  "Border change", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(48 82  "Government defected", place(e) size(small) orientation(horizontal) col(black) )  
		ttext(48 93  "Government changed", place(e) size(small) orientation(horizontal) col(black) )  
	;
	# delimit cr
	
	gr_edit .plotregion1.bars[5].draw_view.setstyle, style(no)
	gr_edit .plotregion1.bars[1].draw_view.setstyle, style(no)
	gr_edit .scaleaxis.reset_rule 0 70 1 , tickset(major) ruletype(range) 
	gr_edit .scaleaxis.style.editstyle majorstyle(tickstyle(show_labels(no))) editcopy
	gr_edit .scaleaxis.style.editstyle majorstyle(tickstyle(show_ticks(no))) editcopy
	
	graph export  "figures/Figure7.pdf", replace
